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Abstract 

Precise knowledge of X-ray diffraction profile shape is crucial in the investigation 
of the properties of matter in crystals powder. Line-broadening analysis is a pre- 
processing step in most of the full powder pattern fitting softwares. Final result of 
line-broadening analysis strongly depends on preliminary three steps: Noise filtering, 
removal of background signal and peak fitting. In this work a new model indepen- 
dent procedure for two of the aforementioned steps (background suppression and 
peak fitting) is presented. The former is dealt with by using morphological mathe- 
matics, while the latter relies on the Hankel Lanczos Singular Value Decomposition 
technique. Real X-ray powder diffraction (XRPD) intensity profiles of Ceria sam- 
ples are used to test the performance of the proposed procedure. Results show the 
robustness of this approach and its capability of efficiently improving the disentan- 
gling of instrumental broadening. These features make the proposed approach an 
interesting and user-friendly tool for the pre-processing of XRPD data. 
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1 Introduction 



X-ray powder diffraction (XRPD) technique is nowadays a well known tool 
to study crystalline properties, which provide important information for ap- 
plications in fields such as nanotechnology [1,2]. All the applications benefit 
from a reliable pre-processing aimed at enhancing the quality of XRPD data. 
The pre-processing procedure consists of four steps: Denoising, background 
suppression, peak fitting and signal deblurring, also known as line-broadening 
(see Fig. 1). Many techniques can be apphed to remove noise from XRPD 
data (sec ref. [3] and references therein). The step of background suppression 
is needed to emphasize the peak features of the sample. Traditional techniques 
are: Young's polynomials, Chebyshev approximation and linear interpolation 
[4]. Peak fitting is a very challenging step to extract information about the 
properties of polycrystalline powder [5]. The final stage in XRPD data pre- 
processing is the signal deblurring aiming at disentangling the instrumental 
line-broadening out of data. On this topic we spot a recent review in ref. [6] . In 
this paper we propose a new approach for two of the above pre-processing steps 
of XRPD data. The signal background is removed by means of a procedure 
based on morphological mathematics. The peak shape profile fitting is carried 
out by using a subspace-based parameter estimation method called Hankel 
Lanczos Singular Value Decomposition (HLSVD) technique [7]. The main ad- 
vantage of this approach is twofold: While both background suppression and 
peak shape profile fitting are model independent, the signal deblurring pro- 
cedure, employed thereafter, crucially benefits from the model independence 
of the previous two steps. Indeed, our signal deblurring steps over the ma- 
jor drawback of most of the XRPD analysis tools available in the literature: 
The peak overlapping problem, namely the difficulty in singling the peak out 
of the full XRPD profile. This problem is mainly due to the model-induced 
bias in background/peak profile reconstruction. As to the denoising we use 
a wavelet based filter, a popular tool available in several software packages. 
Real XRPD intensity profiles of Ceria samples [6] are used to test the per- 
formances. Four different raw datasets were used in pairs [8]. For each pair, 
one dataset was collected on the annealed Ceria specimen (representing the 
instrumental broadening) and the other was collected on the broadened sam- 
ple. The selected pairs are those measured at the University of Birmingham 
(a high resolution X-ray laboratory) and at the National Synchrotron Light 
Source (NSLS X3B1). 
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The structure of the paper is as follows. Section 2 is devoted to the description 
of the proposed method and of the results. Conclusions are drawn in Section 
3. 



2 The method 

Two different raw datasets were downloaded for the Ceria sample: The in- 
strumental standard representing the instrumental broadening and the Ceria 
XRPD pattern of the broadened sample. 

Our strategy in disentangling the profile broadening out of the experimental 
sample relies on a three step procedure which is sketched in the sequel. 

2.1 Denoising 

The noise was removed by applying wavelet transforms to the full XRPD 
spectrum and subtracted prior to the background suppression. 

Among the many applications of wavelets, signal denoising has been deeply 
investigated and the wavelet filter can be considered as the state of art on this 
subject. The discrete wavelet transform (DWT) is a linear operator which 
modifies the data vector in a similar way as the discrete Fourier transform 
(DFT). In both cases, the transform, given hy a, N x N matrix acting on 
the input A^"- vector data, is invertible [9]. The matrix entries are combina- 
tions of basis functions (the familiar sines and cosines in the case of DFT). 
Tables 1-2 summarize the main properties of some wavelet bases. An interest- 
ing property of wavelet basis is the localization in both space and frequency 
domain. This means that they have a finite support or a decay in both do- 
mains. Regularity is another important property of wavelet basis. Regularity 
r means that the r*^ derivative exists alm,ost everywhere (see Table 1). In 
the case of Daubechies wavelet, regularity depends on the order N (see Table 
2). A wavelet is defined by particular set of numbers {cfcj^^Q 27v-i' called 
wavelet filter coefficients. They are determined by imposing the constraints of 

2JV-1 

N vanishing moments, ^ {—l)^k"^Ck = (m = 0, . . . , 2 A" — 1), see Table 

fc=0 

27V-1 

1, and orthogonality ^ CkCk+2m = 25o,m- The coefficients characterize a 

k=0 

low-pass filter while the coefficients bk = {—l)^C2N-i-k results in a high-pass 
filter. Coefficients Ck give the entries oi a N x N matrix, which iteratively 
applies to the A^-data vector thus resulting in the A'-vector of detail coeffi- 
cients. The whole procedure described above is the wavelet transform (sec [9] 
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for further details). Generally speaking, the denoising procedure involves three 
steps. The basic version of the procedure is the following: 

• Calculate the wavelet transform of XRPD profile and sort the components 
of the output vector by increasing frequency. This shall result in A'"-vector 
containing the XRPD profile average coefficient and a set of detail coeffi- 
cients. 

• Noise thresholding, calculated on the highest frequency detail coefficient of 
the wavelet spectrum. 

• Signal reconstruction by using the average coefficient and thresholded detail 
coefficients. 

In this paper we choose the Daubechies wavelet basis with N — 2. We address 
the reader to Daubechies [10] for further details. 



2.2 Background suppression 



The background was determined by means of morphological transforms for 
the full XRPD spectrum and subtracted prior to the signal deblurring. 

The morphological mathematics is based on the language of set theory. Con- 
sidered a discrete binary image X e and a structuring element <S e 3?^, the 
four basic morphological mathematical operations on X by 5 are: 

Dilation: J© 5 = Uses^s 

Erosion: TqS = Hg^s^-s 

Opening : IOS^(IeS)®S 

Closing: I»S ^ {I ® S) Q S 



where Ig denotes the translation of X by s, namely X^ = {a; + s\x G X}. The 
value of each pixel in the output image is based on a comparison of the cor- 
responding pixel in the input image with its neighbours, whose number and 
location is given by the structuring element. Generally, dilation expands im- 
age objects, whereas erosion shrinks them. In practice, dilation and erosion 
are employed in pairs. Opening is the erosion of an image followed by the 
dilation of the eroded image, and closing is the dilation of an image followed 
by the erosion of the dilated image. Opening eliminates sharp peaks smaller 
then the structuring element while the closing fills in the small holes and gaps. 
The binary morphological operations of dilation, erosion, opening and closing 
can be extended to grey-scale images. Let Sj and £s be the domains of the 
gray-scale image X and the grey-scale structuring element S respectively. The 
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grey-scale dilation and erosion can be computed by 

Dilation : (X © 5) {x, y) = max{T{x — m,y — n) + S{m, n)} 

(2) 

Erosion : {T Q S) {x, y) — min {T{x — m, y — n) — S{m, n)} 

where {x — m,y — n) e Ex and (m, n) e Es- For such images, the minimum 
and maximum values are computed within neighbourhood represented by the 
structuring element (see [11] for details). 

In our background suppression procedure the XRPD pattern is reshaped and 
padded into a 2-D image. A disk with a radius of three pixels is used as struc- 
turing element both for erosion and for dilation. As to the erosion (dilation), 
pixels beyond the image border are assigned the maximum (minimum) value 
afforded by the data type. The morphological opening removes small objects 
from the image while preserving the shape and size of larger objects in the 
image. The overall result is a peak smearing effect while the background in- 
tensity remains unalterated. Restoring the original 1-D pattern provides the 
XRPD spectrum background. Figure 2 summarizes the whole procedure. We 
compared our findings to the traditional interpolation method and we found 
a satisfactory agreement (the percentage difference between the background 
computed by traditional techniques and our finding is below 3 %). Up to our 
knowledge, this technique has never been applied to XRPD spectrum back- 
ground suppression and it provides a reliable and user independent estimate 
of it. 



2.3 Peak fitting 



The main problem in analysing an XRPD spectrum is the peak search, since 
the exact d position is crucial in the extraction of the relevant microstructural 
information. Had the peak well defined, its shape would be straighforwardly 
achieved (for instance by a high resolution interpolation/fit by means of a 
model) . Unfortunately the data resolution is rarely high enough to reach the 
goal of a well profiled peak. In that respect several methods have been devised 
so far to accomplish the peak fitting by means of gaussian, lorentzian, voigt, 
pseudo-voigt, Pearson VII and other models [5]. The main drawback of the 
aforementioned methods cited above is the dependence of results on the model 
used in the fit procedure itself and the poor description of the real peak profile 
shape (for instance asymmetry). Here we use HLSVD method to the purpose 
[7]. The main advantage of this method is the fiexibility since the number of 
parameters is not fixed and it can be chosen to achieve a more satisfactory 
agreement between the model and the real peak profile shape. The HLSVD 
method works as follows. Let us model the XRPD intensity samples col- 
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lected at angles n — 0, . . . N — 1 as the sum of K exponentially damped 
complex sinusoids 



K 



In^^Clk exp(-4 -dn) C0S[2'K fk-dn + ^k] , 
k=l 



(3) 



Ofc is the amplitude, (pk the phase, dk the damping factor and fk the fre- 
quency of the k*^ sinusoid, k — 1,...,K, with K the number of damped 
sinusoids. The N data points defined in (3) are arranged into a Hankel ma- 
trix H Hfjj^ = Im+u m = 0, . . . , M - 1, I = 0, . . . , L - 1 , with 
L+M = N+1 {M c:iLc^ N/2). The SVD of the Hankcl matrix is computed as 
Hlxm = Ulxl^lxmV^j^m, where S = diag{Ai, A2, . . . , A,.}, Ai > A2 > . . . A,., 
r = min(L, M), U and V are orthogonal matrices and the superscript H de- 
notes the Hermitian conjugate. The Lanczos bidiagonalization algorithm with 
partial reorthogonahzation is used to compute SVD. This algorithm, based 
on FFT, computes the two matrix-vector products which are performed at 
each step of the Lanczos procedure in 0{{L + M)log2{L + M)) rather than in 
0{LM). In order to obtain the "signal" subspace, the matrix H is truncated 
to a matrix Hk of rank K Hk — Uk^kV^ , where Uk, Vk, and T^k are de- 
fined by taking the first K columns of U and V, and the K x K upper-left 
matrix of E, respectively. As subsequent step, the least-squares solution of the 
following over-determined set of equations is computed vj^°^^E^ ~ yjfiottom) ^ 
where yj^"^^""^^ ^^d vj^"^'' are derived from Vk by deleting its first and last 
row, respectively. The K eigenvalues Zk of matrix E are used to estimate the 
frequencies fk and damping factors dk of the model damped sinusoids from 
the relationship 



with k = 1, . . . ,K. Values so obtained are inserted into the model equation 
(3) which yields the set of equations 



with n = 0, . . . , A^ — 1. The least-squares solution of (5) provides the amplitude 
dk and phase (fk estimates of the model sinusoids which are used in the next 
step. 




(4) 



K 



In-^dk exp(-4 {}n) COs[27r/fc1?n + <Pk] , 
k=l 



(5) 
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2-4 Deblurring 



The XRPD pattern has to be corrected for the instrumental broadening. Sev- 
eral methods have been devised so far to deal with this problem. Among them 
we quote the Stokes method [12] and the Bayesian approach [13,14]. The main 
drawbacks of these methods stem from the difficulty in evaluating the back- 
ground level, mainly due to peak overlapping. 

The XRPD pattern plugged in the deblurring algorithm is noise-background 
free since it has been already pre-processed by the wavelets filter and the 
morphological operator. 

The technique proposed in this paper is a modified version of the one pre- 
sented in ref . [6] . A blurred or degraded XRPD pattern can be approximately 
described by a Voltcrra equation g — T-C(E) f + n, where g is the blurred XRPD 
pattern and Ti is the distortion operator due to several causes, namely the 
point spread function (PSF), and n is an additive noise, introduced in the 
XRPD acquisition, that corrupts the signal. Strictly speaking in an XRPD 
experimental setup we deal with a poissonian noise which is a multiplicative 
noise. However the Poisson distribution function resembles the Gauss one pro- 
vided a sufficiently large statistics in photons counting. 

As to the deblurring procedure we implement the damped Lucy-Richardson 
algorithm. This function performs multiple iterations, using optimization tech- 
niques and Poisson statistics. In our approach the PSF is the raw dataset 
downloaded for the Ceria sample - the instrumental standard - resembling the 
instrument profile [6]. The algorithm maximizes the likelihood that the result- 
ing image, when convolved with the PSF, is an instance of the blurred image, 
assuming Poisson noise statistics. This function can be effective when the PSF 
is known but the knowledge about the additive noise in the image is poor. 
The Lucy-Richardson algorithm introduces several adaptations to the origi- 
nal maximum likelihood algorithm that addresses complex image restoration 
tasks. By using these adaptations, the effect of noise amplification on image 
restoration can be reduced, nonuniform image quality can be accounted for 
(e.g., bad pixels, flat-field variation) and the restored image resolution can be 
improved by subsampling. 

Due to the denoising/background suppression, the original Volterra equation 
is readily simplified: g = Ti® f , where Ti and / have been already defined. 
Was the inverse Ti^^ explicitly known, we would solve the former equation 
at a glance. Unfortunately this is not the case and the solution has to be 
approximated as follows, g, Ti. and / are positive and the PSF cannot change 
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the norm, i.e. \ \g\ \ — ||/||. Thus: 

E/^ = E^^ ^ E/^ - T.9r%^ = T.fr i9®n , (6) 



where Ofe = ,^ p^ ] the solution can be found by an iterative procedure: 

/)fc 

j-n+i _ rn where the initial ffuess XRPD spectrum is uniform. As al- 

ready stressed in ref . [6] , the main drawbacks in applying such an algorithm to 
the single peak deconvolution are the noise amplification and the peak fitting 
bias. Noise amplification is dramatically reduced by both the denoising pro- 
cedure and the small (some five) number of iterations used in the algorithm. 
Moreover the damp in the algorithm specifics the threshold level for the de- 
viation of the resulting image from the original image, below which damping 
occurs. For pixels that deviate in the vicinity of their original values, iterations 
are suppressed. 

As to the peak fitting bias, unlike the Balzar approach, our procedure uses 
the instrumental standard pattern (with no overlapping) as the PSF to decon- 
volve the full XRPD pattern and then we extract the deconvoluted/deblurred 
XRPD pattern in the same range of the PSF used for the deconvolution it- 
self. The rationale of this choice relies on the fact that while the PSF peaks 
have no overlap, this is not the case for the broadened sample peaks and, 
thus, the peak ranges can be defined starting on the annealed sample rather 
than the broadened one. Moreover the discrete Fourier transform (DFT), used 
by the deblurring functions, assumes that the frequency pattern of an image 
is periodic. This assumption creates a high-frequency drop-off at the edges 
of an overlapping peaks cluster [15]. This high-frequency drop-off can create 
an effect called boundary related ringing in deblurred images, that is a sys- 
tematic error affecting any further investigation on the physical meaning of 
the deconvolved spectrum. To reduce ringing, our full pattern deconvolution, 
as described above, resembles an edgetaper function. It removes the high- 
frequency drop-off at the edge of an image by blurring the entire image and 
then replacing the center pixels of the blurred image with the original one. In 
this way, the edges of the image taper off to a lower frequency. 



3 Conclusions 



In this paper wc presented a new approach for background removal and peak 
fitting of XRPD profiles. Such operations are crucial in the line-broadening 
analysis of X-ray diffraction profiles, an important pre-processing step in the 
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investigation of the crystal powder samples by means of XRPD data. Back- 
groud suppression relies on the use of morphological mathematics while peak 
fitting is carried out by means of HLSVD technique. In order to enhance the 
signal-to-noise ratio of XRPD profiles a wavelet based filter is prehminarly 
applied to XRPD data. The output of the proposed precedure is supphed to 
a damped Lucy-Richardson algorithm for deblurring. The main advantage of 
this approach is twofold: Background suppression and peak shape profile fit- 
ting are model independent. Real XRPD intensity profiles of Ceria samples 
are used to test performances. Results show that the background estimate is 
in agreement with that computed by traditional interpolation methods with 
a percentage difference below 3%. Further, the output XRPD profile has nar- 
rower peaks, located at the same position and with the same shape as the 
original one. It is worth noting that once the deblurred, noise-background free 
XRPD spectrum is convoluted back to the PSF and added to the removed 
noise and background signals, the resulting XRPD pattern resembles the orig- 
inal one. 
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Table 1 

Properties of some wavelet bases. 



wavelet 


t-localization 


f-localization 


# zero moments 


r 


Haar 


[0,1] 


1/f 


1 





Sine 


1/t 


[0,1] 


oo 


oo 


Daulxx'hies (X) 


[0,2X-1] 


1/r 


X 


n(X) 



Table 2 

Regularity r = a{N) of a Daubechies wavelet basis of order N. 



N 


2 


3 


4 


5 


6 


N 


a(N) 


0.500 


0.915 


1.275 


1.596 


1.888 


0.2075 N 
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Fig. 1. Block diagram of the pre-processing procedure. The proposed approach refers 
to the background suppression and peak fitting steps (shadowed boxes). 
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Fig. 2. Background evalutation procedure. Left top: original XRPD pattern. Right 
top: XRPD pattern after two-dimensional reshaping. Right bottom: two-dimensional 
reshaped XRPD pattern after morphological operations. Left bottom: background 
after one-dimensional reshaping. 
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Fig. 3. Plots. Left top: original XRPD pattern. Right top: XRPD pattern after 
denoising and background suppression. Left bottom: noise-background free XRPD 
pattern after deblurring. Right bottom: residue between the final XRPD pattern 
re-convoluted to the PSF together with noise and background and the original 
XRPD pattern. 
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